/*
 *  elasticlog.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde Martin on 04/05/11.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "ogden.h"

Ogden::Ogden()
:	HyperelasticStretches(),
	N(3),
	b(0.0)
{}

// ***

Ogden::Ogden(ifstream &in, const string& na, const string & t)
:	HyperelasticStretches(na,t),
 	N(3),
	b(0.0)
{
	int i, j;
	string keyword;
	string line;
	vector<string>	alphaValues;
	stringstream stnum;
	double value;
	
	while (in >> keyword && keyword.compare("ENDMATERIAL") )
	{
		if(!keyword.compare("VOL"))
		{
			in >> typeVol;
			if(!typeVol.compare("mixed")) in >> b;
		}
		//in >> keyword;
		if(!keyword.compare("ALPHA"))
		{
			getline(in,line);
			istringstream iss(line);
			while (iss >> keyword )
			{
				//cout << "keyword is " << keyword << endl;
				alphaValues.push_back(keyword);
			}
			N = alphaValues.size();
			for(i=0; i<N; ++i)	// Starts from 1 because the first element is the keyword V0
			{
				stnum << alphaValues[i];
				stnum >> value;
				vAlpha.push_back(value);
				stnum.clear();
			}
		}
		//in >> keyword;
		if(!keyword.compare("PARTIALMU"))
		{
			for (i=0; i<N; ++i)
			{
				in >> value;
				vMu.push_back(value);
			}
		}
		if(!keyword.compare("DENSITY"))	in >> density;
		if(!keyword.compare("K"))	in >> K;
	}

	mu = 0.0;
	for(i=0; i<N; ++i)	mu += vAlpha[i] * vMu[i] / 2.0;
	
	E = 9.0 * K * mu / (3.0 * K + mu);
	
	lambda = K - 2.0 / 3.0 * mu; 
	
	nu = (3.0 * K - 2.0 * mu) / (2.0 * (3.0 * K + mu));
	
	hypertype = "decoupled";

	Svol.setZero(6);
	Siso.setZero(6);
	Ctvol.setZero(6,6);
	Ctiso.setZero(6,6);


	//cout << "beta= " << b << endl;

	//cout << "E= " << E << endl << "nu= " << nu << endl << "mu= " << mu << endl << "rho= " << density << endl;
	//printSTLVector(vAlpha);
	//printSTLVector(vMu);
}

// ***

void Ogden::strainEnergy()
{
	int i;
	if(J > 0.0)
	{
		Wiso = 0.0;
		for(i=0; i<N; ++i)
		{
			Wiso += vMu[i] / vAlpha[i] * (pow(lambda1m, vAlpha[i]) - 1.0) +
			        vMu[i] / vAlpha[i] * (pow(lambda2m, vAlpha[i]) - 1.0) +
			        vMu[i] / vAlpha[i] * (pow(lambda3m, vAlpha[i]) - 1.0);
		}
		
		if(!typeVol.compare("cuad"))	Wvol = K / 2.0 * (J - 1.0) * (J - 1.0);
		if(!typeVol.compare("log"))		Wvol = K / 2.0 * log(J) * log(J);
		if(!typeVol.compare("mixed"))	Wvol = K / (b * b) * (1.0 / pow(J,b) - 1.0 + b * log(J));
		
		W = Wiso + Wvol;
	}
	else
	{
		cout << "*ERROR: Negative volume, material " << name << " has collapsed \n";
		Wvol = -0.5;
		Wiso = -0.5;
	}
}

// ***

void Ogden::DstrainEnergy()
{
	int i;
	if(J > 0.0) 
	{ 
		Wiso1 = 0.0;
		Wiso2 = 0.0;
		Wiso3 = 0.0;
		
		for(i=0; i<N; ++i)
		{
			Wiso1 += vMu[i] * pow(lambda1m, vAlpha[i] - 1.0);
			Wiso2 += vMu[i] * pow(lambda2m, vAlpha[i] - 1.0);
			Wiso3 += vMu[i] * pow(lambda3m, vAlpha[i] - 1.0);
		}

		if(!typeVol.compare("cuad"))	dWvol = K * (J - 1.0);
		if(!typeVol.compare("log"))		dWvol = K * log(J) / J;
		if(!typeVol.compare("mixed"))	dWvol = K / (b * J) * (1.0 - 1.0 / pow(J,b));
	}
	else
	{
		dWvol = - 1.0;
		Wiso1 = - 1.0;
		Wiso2 = - 1.0;
		Wiso3 = - 1.0;
	}
}

// ***

void Ogden::DDstrainEnergy()
{
	int i;
	if(J > 0.0)
	{
		Wiso11 = 0.0;
		Wiso12 = 0.0;
		Wiso13 = 0.0;
		Wiso22 = 0.0;
		Wiso23 = 0.0;
		Wiso33 = 0.0;
		
		for(i=0; i<N; ++i)
		{
			Wiso11 += vMu[i] * (vAlpha[i] - 1.0) * pow(lambda1m, vAlpha[i] - 2.0);
			Wiso22 += vMu[i] * (vAlpha[i] - 1.0) * pow(lambda2m, vAlpha[i] - 2.0);
			Wiso33 += vMu[i] * (vAlpha[i] - 1.0) * pow(lambda3m, vAlpha[i] - 2.0);
		}
		
		if(!typeVol.compare("cuad"))	ddWvol = K; 
		if(!typeVol.compare("log"))		ddWvol = K * (1.0 - log(J)) / (J * J);
		if(!typeVol.compare("mixed"))	ddWvol = K / (b * J * J) * ( (b + 1.0) / pow(J,b) - 1.0);
	}

	else
	{
		ddWvol = - 1.0;
		Wiso11 = - 1.0;
		Wiso12 = - 1.0;
		Wiso13 = - 1.0;
		Wiso22 = - 1.0;
		Wiso23 = - 1.0;
		Wiso33 = - 1.0;
	}
}

// ***

void Ogden::DDDstrainEnergy()
{
	int i;
	if(J > 0.0)
	{
		Wiso11 = 0.0;
		Wiso12 = 0.0;
		Wiso13 = 0.0;
		Wiso22 = 0.0;
		Wiso23 = 0.0;
		Wiso33 = 0.0;
		
		for(i=0; i<N; ++i)
		{
			Wiso11 += vMu[i] * (vAlpha[i] - 1.0) * pow(lambda1m, vAlpha[i] - 2.0);
			Wiso22 += vMu[i] * (vAlpha[i] - 1.0) * pow(lambda2m, vAlpha[i] - 2.0);
			Wiso33 += vMu[i] * (vAlpha[i] - 1.0) * pow(lambda3m, vAlpha[i] - 2.0);
		}
		
		if(!typeVol.compare("cuad"))	ddWvol = K;
		if(!typeVol.compare("log"))		ddWvol = K * (1.0 - log(J)) / (J * J);
		if(!typeVol.compare("mixed"))	ddWvol = K / (b * J * J) * ( (b + 1.0) / pow(J,b) - 1.0);
	}
    
	else
	{
		ddWvol = - 1.0;
		Wiso11 = - 1.0;
		Wiso12 = - 1.0;
		Wiso13 = - 1.0;
		Wiso22 = - 1.0;
		Wiso23 = - 1.0;
		Wiso33 = - 1.0;
	}
}

// ***

//double Ogden::getE()
//{	return E;}
//
//// ***
//
//double Ogden::getNu()
//{	return nu;}
//
//// ***
//
//double Ogden::getLambda()
//{	return lambda;}
//
//// ***
//
//double Ogden::getMu()
//{	return mu;}
//
//// ***
//
//double Ogden::getDensity()
//{	return density;}

// ***

vector < double > Ogden::getvMu()
{	return vMu;}

// ***

vector < double > Ogden::getvAlpha()
{	return vAlpha;}

// ***

//double Ogden::getBulkmodulus()
//{	return K;}
